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In the relativistic and the nonrelativistic theoretical treatment of moderate and high-power laser- 
matter interaction, the generalized Bessel function occurs naturally when a Schrodinger-Volkov 
and Dirac-Volkov solution is expanded into plane waves. For the evaluation of cross sections of 
quantum electrodynamic processes in a linearly polarized laser field, it is often necessary to evaluate 
large arrays of generalized Bessel functions, of arbitrary index but with fixed arguments. We show 
that the generalized Bessel function can be evaluated, in a numerically stable way, by utilizing a 
recurrence relation and a normalization condition only, without having to compute any initial value. 
We demonstrate the utility of the method by illustrating the quantum-classical correspondence of 
the Dirac-Volkov solutions via numerical calculations. 
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I. INTRODUCTION 



form exp [ix s'md — iy sin(20)] as follows. 



The Volkov solution [1] is the exact solution of the 
Dirac equation in the presence of a classical plane-wave 
laser field of arbitrary polarization. In order to evaluate 
cross sections by quantum electrodynamic perturbation 
theory, it is crucial to decompose the Volkov solutions 
into plane waves, in order to be able to do the time and 
space integrations over the whole Minkowski space-time. 
If the laser field is linearly polarized, one naturally en- 
counters the generalized Bessel functions as coefficients 
in the plane-wave (Fourier) decomposition of the wave 
function, both for the Dirac-Volkov equation as well as 
for the laser-dressed Klein-Gordon solutions, and even 
for Schrodinger-Volkov states (see also Sec. W} . 

The wide use of the generalized Bessel function in the- 
oretical laser physics is thus due to the fact that different 
physical quantities, such as scattering cross sections and 
electron-positron pair production rates, can be expressed 
analytically in terms of infinite sums over generalized 
Bessel functions which we here denote by the symbol 
Jn(x,y). The generalized Bessel function Jn(x,y) is a 
generalization of the ordinary Bessel function Jn{x) and 
characteristic of the interaction of matter with a linearly 
polarized laser field; it depends on two arguments x and 
y, and one index n. Here, we use it in the convention 



1 r 

Jn{x,y) = — / cxp[— in0 + ix sin(0) 
27r 



(1) 



ly sin(26l)] d6l. 



where n is an integer, and x and y are real numbers. 
Jn(x,y) is real valued. The generalized Bessel functions 
provide a Fourier decomposition for expressions of the 



exp [ia; sin0 — iy sin(26')] = J„(a;, y) exp (i n 6*) . 

n— — oo 

(2) 

In practical applications, the angle 6 often has the physi- 
cal interpretation of a phase of a laser wave, 9 = ut—k-x, 
where uj is the angular laser photon frequency, and k is 
the laser wave vector. By contrast, the well-known ordi- 
nary Bessel functions are defined as 

1 r 

Jn[x) — — / exp [— in6' -f ix sin(6')] d6' , (3) 
and they have the fundamental property 

oo 

exp {ix sin 9)— ^ Jn{x) exp{in9) . (4) 



n— — oo 
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The generalized Bessel function was first introduced by 
Reiss in the context of electron-positron pair creation , 
followed by work of Nikishov and Ritus and Brown 
and Kibble Further examples of work utilizing 

Jn{x,y) in the relativistic domain include pairproduc- 
tion by a Coulomb field and a laser field d, H, Q , laser- 
assisted bremsstrahlung [^, muon-antimuon cre- 
ation J^], undulator radiation [l^, and scattering 
problems, both classical [O] , and quantum mechani- 
cal [H, [l3| . A fast and reliable numerical evaluation of 
Jn{x,y) would also speed up calculation of wave packet 
evolution in laser fields EBl- In nonrelativistic cal- 
culations, the generalized Bessel function has been em- 
ployed mainly for strong-field ionizatio n fiol. |20| . 
but also for high-harmonic generation |23l. |24{ . 

On the mathematical side, a thorough study of Jn{x, y) 
has been initiated in a series of papers [2^, [2^, [27j . 
and even further generalizations of the Bessel function 
to multiple arguments and indices have been considered 
[28l . 12^ . IsO . |31| . On the numerical side, relatively little 
work has been performed. Asymptotic approximations 



2 



have been found for specific regimes and a uniform 

asymptotic expansion of Jn{x,y) for large arguments by 
saddle-point integration is developed in Ref. [s^- For 
some of the applications described above, in particu- 
lar when evaluating second-order laser-assisted quantum 
electrodynamic processes [33l |. a crucial requirement is 
to evaluate large sets of generalized Bessel functions, at 
fixed arguments x and y, for all indices n for which 
the generalized Bessel functions acquire values which 
are numerically different from zero (as we shall see, for 
|n| » the generahzed Bessel functions decay ex- 

ponentially with n). 

It is clear that recursions in the index n would greatly 
help in evaluating large sets of Bessel functions. For or- 
dinary Bessel functions, an efficient recursive numerical 
algorithm is known, and it is commonly referred to as 
Miller's algorithm [s^, [sHi- However, a generalization of 
this algorithm for generalized Bessel functions has been 
lacking. The purpose of this paper is to provide such 
a recursive numerical algorithm: We show, using ideas 
from [H, [13, [H, di'l , that a stable recurrence algorithm 
can indeed be established, despite the more complex re- 
currence relation satisfied by J„(x,?/), as compared to 
the ordinary Bessel function Jn{x). The reduction of 
five-term recursions to four- and three-term recursions 
proves to be crucial in establishing a numerically stable 
scheme. 

The computational problem we consider is the follow- 
ing: to evaluate 



Jn(x,y) : X fixed, y fixed, 
where rimin <n< n^ax , 



(5) 



by recursion in n. Our approach is numerically stable, 
and while all algorithms described here have been im- 
plemented in quadruple precision (roughly 32 decimals), 
we note that the numerical accuracy of our approach can 
easily be increased at a small computational cost. 

Our paper is organized as follows. In Sec. |TT1 we recall 
some well-known basic properties of J„(x,?/), together 
with some properties of the solutions complementary to 
Jn{x,y), which fulfill the same recursion relations (in n) 
as the generalized Bessel functions but have a different 
asymptotic behavior for large \n\ as compared to Jn{x, y). 
After a review of the Miller algorithm for the ordinary 
Bessel function, we present a recursive Miller-type algo- 
rithm for generalized Bessel functions in Sec. IIIIl and 
show that it is numerically stable. In Sec. IIVI we nu- 
merically study the accuracy which can be obtained, and 
compare the method presented here with other available 
methods. We also complement the discussion by con- 
sidering in Sec. IVl illustrative applications of the numer- 
ical algorithm for Dirac-Volkov solutions in particular 
parameter regions, together with a physical derivation of 
the recurrence relation satisfied by the generalized Bessel 
function. Section IVll is reserved for the conclusions. 



II. BASIC PROPERTIES OF THE 
GENERALIZED BESSEL FUNCTION 

A. Orientation 

Because the definition ^ provides us with a conve- 
nient integral representation of the generalized Bessel 
function, all properties of Jn(x,y) needed for the follow- 
ing sections of this article can in principle be derived from 
this representation alone fs','!^. E.g., shifting 9 —0—tt 
and 9 9 + Tr, respectively, in ([1]) gives two symmetries. 



Jn{x, -y) 
Jn{-x,y) 



-l)"J„„(x,y), 
'irJn{x,y), 



(6) 



from which J_„(x,?/) = Jn{—x,—y) follows. We recall 
the corresponding properties of the ordinary Bessel func- 
tion. 



Jn{x) = i-irJni-x) = (-1)" J_„(x) , 



(7) 



Due to the symmetries ©, we can consider in the fol- 
lowing only the case of positive x and y without loss of 
generality, provided we allow n to take arbitrary posi- 
tive and negative integer values. Our sign convention 
for the y sin 20-term in the argument of the exponential 
in Eq. ([T]) agrees with but differs from the one used 
in [1^. As is evident from inspection of Eqs. ^ and ([3]), 
Jn{x, y) can be expressed as an ordinary Bessel function 
if one of its arguments vanishes. 



_ / J-n/2[y) if n even 







if n odd. 



J„(x, 0) = J„(x), Jn(0,y) 

I VJ Li. I L \J\JL 

(8) 

By inserting the expansion of the ordinary Bessel func- 
tion Yl'^=-tx, Jn{x) exp{in9) = exp(ix sui9) into Eq. ([1]), 
we see that Jn(x,y) can be expressed as an infinite sum 
of products of ordinary Bessel function, 

C30 

Jn{x,y)^ ^ J2s+n{x)Js{y)- (9) 

S— — 00 

There are also the fohowmg sum rules, 

oo oo 

Jn{x,y)= Jl{^,y) = ^, (10) 



n— — oo 



n— — oo 



which can be derived by considering the case 9 = Q in 
Eq. (O [for Yl'^=-oGJri{x,y) = 1], and by considering 
Eq. (I2|) multiplied with its complex conjugate, and inte- 
grating over one period [for Yl^=-Qa •^ni^^v) ^ The 
relation (jlOp is important for a recursive algorithm, be- 
cause it provides a normalization for an array of gener- 
alized Bessel function computed according to the recur- 
rence relation 



2nJn{x, y) =x [J„+i(a;, y) + Jn-iix, y)] 
-2y[Jn+2(x,y) + Jn-2{x,y)] , 



(11) 
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which connects generahzed Bessel functions of the same 
arguments but different index n. Equation (jlip can be 
derived by partial integration of Eq. ([1]). Interestingly, 
Eq. pT|) together with the normalization condition (fTU|) 
can be taken as an alternative definition for Jn{x, y), from 
which the integral representation ^ follows. The recur- 
sion pT|) is the basis for the algorithm described below 
in Sec. IIIIl 



B. Saddle point considerations 

A qualitative picture of the behavior of Jn{x,y) as a 
function of n can be obtained by considering the posi- 
tion of the saddle points of the integrand in ^ [13, H^l • 
By definition, a saddle point Os denotes the point where 
the derivative of the argument of the exponential in ([T]) 
vanishes, and therefore satisfies 



e 

By writing Jn{x,y) as 
Jn[x,y) 



± 



1 n 



64 y^ 2 Ay 



(12) 



1 { r 

— Re< / exp[— in0 + i.T sinfS) 

7^ yJo 



ly sin(26l)]d6l 



(13) 



we can consider only saddle points with < Re6's < tt. 
By the properties of the cosine function, the saddle points 
come in conjugate pairs, so that if 9s is a saddle point, 
so is 0*. Furthermore, since cos(27r — dg) = cos^s, the 
saddle points are placed mirror symmetrically around 
lie 9s = TT. Since each of the endpoint contributions at 
9 = and 9 = tt to the integral (fT3| vanish (provided 
the endpoints are not saddle points), an asymptotic ap- 
proximation for Jn(x,y) is provided by the saddle point 
method (the method of steepest descent) ^^\, by sum- 
ming the contributions from the saddle points 9s situated 
on the path of steepest descent. Here, imaginary saddle 
points (i.e., saddle points with lm9s± ^ 0) give expo- 
nentially small contributions to the integral, while real 
saddle points contribute with an oscillating term. Closer 
inspection of Eq. (fT2|) reveals two cases. 

In case 1, with 8t/ > x, there are four different regions, 
which we denote by ai , 5i , ci , di (see Table HI . In region 
fli, where n < —2y—x, we have four distinct saddle points 
solutions 9s±, which are all imaginary, and Jn{x,y) 
is exponentially small. Region bi, where —2y — x < n < 
—2y+x, has two imaginary {9s+, 9*^) and one real saddle 
point 9s-, and Jn{x,y) exhibits an oscillating behavior 
here. For —2y + x < n < 2y + x^/{16y), i.e. in region 
Ci, both saddle points are real, and in the region di, 
n > 2y + x'^/{16y), the two saddle points 9s± are again 
imaginary, which results in very small numerical values 



TABLE I: Saddle-point configurations for the generalized 
Bessel function Jn{x,y) as a function of the arguments x 
and y. A distinct imaginary saddle point is denoted "imag." 
whereas a real saddle point is denoted "real." The different 
regions are illustrated in Fig. [T] 



Case 1: 


> X 




region 


condition 


saddle points 






^ illlfXti. 


fti 
(^1 


'Oqi . — - T T\ <^ - — -'?7( -4- T 


O imao' — l-TP^l 


Cl 


-2y + x <n<2y + x'^/(16y) 


2 real 




n > 2y + x'^/{16y) 


2 imag. 


Case 2: 


8y < X 




region 


condition 


saddle points 


a2 


n < -~2y — X 


4 imag. 


&2 


—2y — X < n < ~2y + x 


2 imag.-|-real 


C2 


n > —2y + X 


2 imag. 



of the generalized Bessel functions. For case 2, 8y < x, 
there are only three regions, as recorded in Table [J The 
two cases coincide if 8y = x. Figure [T] illustrates the two 
different cases. 

In all regions ai, 6i, ci, di, 02, 62, C2, there are, depend- 
ing on the region, up to four saddle points to consider. Of 
these one or two saddle points contribute to the numer- 
ical approximation to Jn{x,y). For large arguments y, 
X and/or a large index n, asymptotic expressions can be 
derived [sol |33 . The general form for the leading-order 
term is (see |40| for a clear exposition of the general the- 
ory of asymptotic expansions of special functions) 




exp (1/(6*8+) -in9s++ ie+) 



exp (i/(6's-) -in9s- + ie-) 



(14) 



F+{n,x,y) 



F-in,x,y), 



where 



f{9) = xsm{9)~ysmi29). 



(15) 



For imaginary saddle points, only the contribution of 
those situated on the path of steepest descent should be 
included , i.e., the integration around the saddle point 
should be carried out along a curve of constant complex 
phase, with 9 = ^ + iv satisfying 



Im(i/(6') -int 



const. 



(16) 
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(a) (b) 

FIG. 1: (Color online.) Illustration of the different saddle point regions of Jn{x,y), for the two qualitatively different cases 
number 1, with 8y > x (here x = y = 10^ was used), and case number 2, with 8y < x {x = lOy = 10'^). In case 1 [panel (a)], the 
transition from region bi, where only one saddle point is real, to ci, where two real saddle points contribute, occurs precisely 
at n = —2y + x = —10^ (see also Table |l]|. The complex oscillating behavior in region Ci can be understood as interference 
between the contributions from the two real saddle points. In case 2 [panel (b)], we have Sy < x, with only three as opposed 
to four qualitatively different regions (see also Table |T]|. 



on that curve. In practice this means for regions with 
imaginary saddle points only, Jn{x,y) is given by the 
contribution from the saddle point with smallest jlm^sj, 
and in regions with both imaginary and real 6s the contri- 
bution from the imaginary saddle point can be neglected. 
However, we will see in the following discussion that all 
saddle points, including those not on the path of steep- 
est descent which would produce an "exponentially large" 
contribution, can be interpreted in terms of complemen- 
tary solutions to the recurrence relation (|lip . The con- 
stant phase e± in Eq. ^1)1 is given by 



■sgn[/"(0.±)] 



(17) 



for real saddle points. For imaginary saddle points, e± 
can be found from the requirement 



dv 

tane± = — 
a/1 



(18) 



with 6 = ^ + iv describing the path of steepest descent 
[see Eq. p6)) ]. For a detailed treatment of the saddle 
point approximation of J„(x,?/) we refer to [s^l, where 
uniform approximations, valid also close to the turn- 
ing points (the borders between the regions described in 
Fig. [1]), and beyond the leading term p4)) . are derived. 
For our purpose, namely to identify the asymptotic be- 
havior of the complementary solutions, the expression 
is sufficient. 



C. Complementary solutions 

The recurrence relation (|lip involves the five general- 
ized Bessel functions of indices n — 2, n— l,n,n+ 1, n + 2. 
In general, an rn-term recursion relation is said to be of 



order m — 1. If we regard the index n as a continuous 
variable, then a recursion relation of order m — 1 corre- 
sponds to a differential equation of order m — 1, which 
has (m— 1) linearly independent solutions. Equation (fTT|) 
consequently has four linearly independent (complemen- 
tary) solutions. The function J„(x,?/) is one of these. 

For the analysis of the recursive algorithm in Sec. IIIII 
below, we should also identify the complementary solu- 
tions to the recurrence relation (|lip . For our purposes, 
it is sufficient to recognize the asymptotic behavior of 
the complementary solutions in the different regions ai- 
di and a2-C2 (see Fig. [1]). It is helpful to observe that 
the recurrence relation pT|) is satisfied asymptotically by 
each term F±{n,x,y) from Eq. (jl4p individually. The 
recurrence relation pT|) is also satisfied, asymptotically, 
by a function obtained by taking in Eq. a saddle 
point that is not on the path of steepest descent, which 
is equivalent to changing the sign of the entire argument 
of the exponential. In addition, for real saddle points and 
in regions with only two imaginary saddle points, the re- 
currence relation is asymptotically satisfied by taking the 
same saddle point but the imaginary part instead of the 
real part in Eq. (fH)) (and thereby changing the phase). 

In regions with four imaginary saddle points (01,02), 
there are thus two solutions that are exponentially in- 
creasing with the index n — > — 00 [the two solutions cor- 
respond to the two saddle points 9s where Re(i/(6's) — 
inOs) > 0], and two further solutions which are exponen- 
tially decreasing [From 9s with Ke{if{9s) — m9s) < 0]. 
In regions with two imaginary and one real saddle point 
(61,62), the four solutions behave as follows. There are 
two oscillatory solutions [these correspond to the real and 
imaginary parts of the term which contains the real sad- 
dle point in Eq. p4p ]. and a third solution which is expo- 
nentially increasing, and a fourth one which is exponen- 
tially decreasing as n — > —00 [the two latter solutions are 
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due to the imaginary saddle points in Eq. ([14])] . The re- 
gion with two real saddle points (ci) has four oscillating 
solutions, as a function of n. Finally, in regions di and C2, 
where we have two distinct imaginary saddle points, we 
have two exponentially increasing (as n oo) solutions 
with different phase, and two exponentially decreasing 
with different phase. Concerning the question of how 
to join the different asymptotic behaviors to form four 
linearly independent solutions, we note that Jn{x,y) is 
the only solution which can decrease in both directions 
n — > ±00, since it represents the unique, normalizable 
physical solution to the wave equation (see subsection 
IV A[) . Furthermore, there must be one solution that in- 
creases exponentially where Jn{x,y) decreases, and that 
exhibits an oscillatory behavior where Jn(a;, y) also oscil- 
lates. The reason is that in either of the limits x ^ 
or y — > 0, we must recover the ordinary Bessel function, 
and the Neumann function as the two solutions to the 
recurrence relation. Having fixed the asymptotic behav- 
ior of two of the solutions, the behavior of the the two 
remaining functions follows. We label the four different 
solutions with J„, Xn, and Zn-, where Jn is the gen- 
eralized Bessel function Jn{x,y) with the arguments x, y 
suppressed. 

Integral representations for the complementary solu- 
tions can be found by employing Laplace's method [4l| . 
details of which will be described elsewhere. The explicit 
expressions can be found in the Appendix. However, as 
noted previously, in this paper we shall need only the 
asymptotic properties of the complementary solutions, 
which can be deduced from p4)) . 

We also observe that the situation for Jn{x,y) de- 
scribed above is directly analogous to that of the or- 
dinary Bessel function Jn{x) and the complementary 
Neumann (also called Weber) function Yn{x) of a sin- 
gle argument. For x ^ n they have the asymptotic 
behavior Jn{x) « Re 1/2/ (ttx) exp(ia; — m/A — im:/2), 
Yn{x) ~ Im •\/2/(7rx) exp(ia; — in/A — imr/2), and for 
x n we have Jn{x) « (ea:)"'(2n)~"/V27rn, Yn{x) ~ 
-2(ea;)-"(2n)"/\/i^. 

According to the above discussion and as illustrated 
in Fig. [21 the functions J„, y„, X„, and Zn have the 
following relative amplitudes in the different regions: 

region ai,a2 : |Z„| > |F„| > |J„| > |X„|, 

region 61,62 : \Zn\ > \Yn\ ^ |J„| > ^^^^ 

region ci : |Z„| |r„| |J„| ~ |X„|, 
region di,C2 : ^ |X„| > |J„| 

In Eq. p9| . we have assumed that all functions have 
the same order of magnitude in the oscillating region. 
This can be accomplished by choosing a suitable constant 
prefactor for the complementary functions Z„, Yn, and 
Xn- Figure [2| shows an example of the four different 
solutions, for case 1 (8y > x). The actual numerical 
computation of the complementary solutions is discussed 
in Sec. HITCl 




n 



FIG. 2: (Color online.) The five-term recurrence relation (|lip 
has four linearly independent solutions. Note the logarithmic 
scale. The values x — y = 10^ were used for the calculation, 
corresponding to case 1 {8y > x). The solutions are labeled 
by Jn [red line, the true generalized Bessel function J„{x,y)], 
Yn (light blue line), Xn (green Hne), and Zn (blue line). The 
numerically obtained solutions X„, Yn, and Zn have been 
shifted vertically by multiplication with an appropriate con- 
stant (of order 10^° for Yn, 10"^" for Zn, and 10"^" for X„) 
for clarity. The separation of the different saddle point re- 
gions is marked with dashed lines. The regions ai,6i,ci and 
di are described in Table |I| 



III. MILLER-TYPE ALGORITHM FOR 
GENERALIZED BESSEL FUNCTIONS 

A. Recursive Miller's algorithm for ordinary Bessel 
functions 

A straightforward implementation of Miller's algo- 
rithm |34i. I35I. |43| can be used for the numerical calcula- 
tion of the ordinary Bessel function Jn{x). We note that 
there are also other ways of numerically evaluating Jn{x), 
which include series expansions .43i] or contour integra- 
tion [3| ■ In the following, we review the simplest form of 
Miller's algorithm, to prepare for the discussion on the 
generalized algorithm. We treat the case of positive n 
and X. For negative n and x, we appeal to the symmetry 
relation ([7]). The properties of Jn{x) used for the algo- 
rithm are the recurrence relation (jlip . with y = Q, which 
automatically reduces (jlip to a three-term relation with 
only two linearly independent solutions. We also use the 
normalization condition Y^^=~oo Jn{x) — 1. 

Viewed as a function of n, Jn{x) exhibits an oscilla- 
tory behavior for n < x, and decreases exponentially for 
n > X. The complementary solution Yn{x), called the 
Neumann function, oscillates for n < x and grows expo- 
nentially for n > X. To calculate an array of Jn{x), for 
< n < N, with > a;, we proceed as follows. We take a 
(sufhciently large) integer M > N , and the initial values 
Cm — 1, cm+1 = 0. We use the recurrence relation fTTI) 
with y — Q to calculate all c„ with indices < n < M 
by downward recursion in n. Now, since the ensemble of 
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the Yn{x) plus the Jn{x) constitute a complete basis set 
of functions satisfying the recurrence relation, the com- 
puted array of the c„ can be decomposed into a linear 
combination, 

c„ = aJ„(x)+/3y„(x), (20) 

where a and /3 are constants, and this decomposition is 
valid for any n. That means that the same decomposition 
must also be valid for the initial index M + 1 from which 
we started the downward recursion, i.e. 

CM+i = = aJM+lix) + PYM+iix). (21) 

From Eq. ^0\j it follows that 

c.=4j.(x)- '^-y(^\^f^ ). (22) 
V Ym+i{x) J 

Provided the starting index A{[ > x is chosen large 
enough, the quantity Jm+i{x) /Ym+i{x) is a small quan- 
tity, due to the exponential character of Jn{x) and Yn{x) 
for index n > a;, so that the computed array c„ is to 
a good approximation proportional to the sought Jn{x). 
Loosely speaking, we can say that we have selected the 
exponentially decreasing function J„ (x) by the downward 
recursion, because the exponentially increasing function 
Yn{x) as \n\ — > oo is suppressed in view of its exponential 
decrease for decreasing |n|. In other words, the error in- 
troduced by the initial conditions decreases rapidly due 
to the rapid decrease of Yn{x) for decreasing n, so that 
effectively only the part proportional to Jn{x) is left. 

Finally, the constant a can be found by imposing the 
normalization condition 

oo 

^c„ = co + 2^C2„ = 1. (23) 

n 71—1 

Here, we have used the symmetries ([7]) in order to elimi- 
nate the terms of odd index from the sum. 

Remarkably, numerical values of J„(x) can be com- 
puted by using only the recurrence relation and the nor- 
malization condition, and not a single initial value is 
needed [e.g., one might otherwise imagine Jo{x) to be 
calculated by a series expansion]. Miller's algorithm has 
subsequently been refined and the error propagation an- 
alyzed b y se veral authors [H, IH, IH, H^l , and also imple- 
mented gi,!!!,!!^. 

B. Recursive algorithm for generalized Bessel 
functions 

In view of the four different solutions pictured in Fig. [21 
it is clear from the discussion in the preceding subsection 
that Jn{x,y) cannot be calculated by naive application 
of the recurrence relation. The general paradigm (see 
Fig. ^ therefore has to change. We first observe that if 
we would start the recursion using the five-term recur- 
rence relation (|lip in the downward direction, starting 



from large positive n, then the solution would eventually 
pick up a component proportional to Z„, which diverges 
for n — oo. Conversely, if we would start the recursion 
using the five-term recurrence relation (jlip in the up- 
ward direction, starting from large negative n, then the 
solution would pick up a component proportional to X„. 
Thus, Eq. PT|) cannot be used directly. 

The solution to this problem is based on rewriting (jlip 
in terms of recurrences with less terms (only three or four 
as opposed to five). By consequence, the reformulated re- 
currence has less linearly independent solutions, and in 
fact it can be shown (see the discussion below) that the 
four-term recurrence, if used in the appropriate directions 
in n, numerically eliminates the most problematic solu- 
tion Zn which would otherwise be admixed to Jn{x,y) 
for n — 5- oo, leading to an algorithm by which it is possi- 
ble to calculate the generalized Bessel function Jn{x,y) 
for n down to the point where we transit from region 
bi to ai in Fig. [2l where the recurrence invariably picks 
up a component from the exponentially growing solution 
y„, and it becomes unstable. However, by using the ad- 
ditional three-term recurrence in suitable directions in 
n, we can numerically eliminate the remaining problem- 
atic solution Yn which would otherwise be admixed to 
Jn{x,y) for n — oo even after the elimination of Z„, 
leading to an algorithm by which it is possible to cal- 
culate the generalized Bessel function J„(x,?/) for n up 
to the point where we transit from region ci to di in 
Fig. [21 where the recurrence invariably picks up a com- 
ponent from the exponentially growing solution X„, and 
it becomes unstable. In the end, we match the results 
of the four-term recursion and the three-term recursion 
at some "matching index" K, situated in region bi or ci, 
normalize the solutions according to Eq. (jTO]), and obtain 
numerical values for Jn{x,y). 

Indeed, in region bi (see Fig. [2]), the wanted solution 
Jn{x,y) satisfies |X„/X„+i| < | J„(a;, y)/J„+i(x, ?;)| < 
\Zn/ Zn+i\^ which means that here application of the re- 
currence relation is unstable in both the upward and 
downward directions with respect to n. By a suitable 
transformations of the recurrence relation, we remove 
one, and then two of the unwanted solutions Yn and 
Zn. With only three (or two) solutions left, we can pro- 
ceed exactly as described in subsection IIII AI to calculate 
Jn{x, y) in a stable way by downward recursion in n. We 
note that the general case of stable numerical solution of 
recurrence relations of arbitrary order has been described 
previously in [H, [13, [H, [H, [5l| , but the application of 
this method to the calculation of the generalized Bessel 
function has not been attempted before, to the authors 
knowledge. 

In the following, we describe the algorithm to com- 
pute an approximation to the array Jn{x, y), Jimin l£ n < 
We let 

2 y -I if 8y>x 

16 y (24) 

~2y + x it 8y < X 
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denote the "cutoflF" indices, beyond which J„{x,y) de- 
creases exponentially in magnitude. In terms of the re- 
gions introduced in TablelH n_ marks the transition from 
region ai to bi for case 1 (or 02 to 62 for case 2), and n_|_ 
is the border between region ci and di for case 1 (be- 
tween 62 and C2 in case 2). Note that we do not assume 
nmin < n- or n,nax > n+, in general r7,,„in and nmax are 
arbitrary (with rimin < f^max)- The usual situation is 
however to require n^i^ < n_ and rimax > n^. With- 
out loss of generality, we assume that both x and y are 
nonzero [otherwise the problem reduces to the calculation 
of ordinary Bessel functions via Eq. ([5])] . 

Central for our algorithm is the transformation of the 
five-term recurrence relation (fTTj) into a four-term and 
three-term recurrence relation. Suppressing the argu- 
ments X and y, we can write the four-term recurrence 



but otherwise arbitrary initial values. A practically 
useful choice is £_li_ — 
four-term formula (|27p and 
the thrcc-tcrm recurrence 



- f3 - 1 for the 
— 1 for 



\2 



2 y + ^nJn+ £.n Jn-1 + Jn-2 



0, 



(25) 



3. Calculate the array /„, with Af_ < n < M+ ac- 
cording to the recurrence formula (j25p . and the ar- 
ray (7„, also with M_ < n < Af_|_ according to the 
recurrence formula (|26p. In both cases, the recur- 
rence is performed in the downward directions for 
n, with arbitrary (but not all zero) starting values 
for /m+, /m++2 and 5m+, 5a/++i- 

4. Choose a "matching index" K , with < K < n+ 
to match the solutions and /„ to each other, 
realizing that gn will be unstable for n > n+, and 
/„ will be unstable for n < n_. Specifically, we 
construct the array 



and the second-order relation 



Jn-1 — . 



(26) 



The coefficients themselves also satisfy recursion rela- 
tions, which are however of first order, namely 



en 
en 

and 



4y^ 

?n-l 
Sn-l 



en - 2(n - 1) 



2?/gLi 

Sn-l 



(27) 



e 

Sn 



2ya 



2 -^n-lSn 



A2 



(28) 



By construction, all sequences ?/„ that solve the orig- 
inal recurrence relation (fTTj) . also solve Eq. and 
Eq. (|26|) . regardless of the initial conditions used to cal- 
culate the coefficients ^^'^'^ and A^'^. The converse does 
not hold: a solution i/„ to the transformed recurrence re- 
lation (P5)) or ([^5]) does not automatically solve Eq. pT|) . 
Rather, this depends on the initial conditions used for 
the coefficients ^'"^'^ i^^ Ki^)- 

The transformation into a four-term and three-term re- 
lation offers a big advantage, as briefly anticipated above. 
We now describe how the algorithm is implemented in 
practice, and postpone the discussion of numerical sta- 
bility imtil subsection nil CI We proceed in five steps. 

1. Select a positive starting index > n+, Timax and 
a negative starting index M_ < n_,7imin, where 
the M's differ from the n's by some "safety mar- 
gin." The dependence of the accuracy obtained on 
the "safety margin" is discussed later, in Sec. IIVI 



2. Calculate the arrays C^'^'^, and A^'^ for M_ < n < 
Af+ 4- 1, employing the recurrence relations (I??)) 
and (|28p in the upward direction of for n. The 



recurrence is started at n = M_ with nonzero j , 



9n 
9K 



K 



if <n< K , 



/„ if K <n< M+ , 



(29) 



where fK,gK^O^s assumed. 



5. The numerical approximation to the generalized 
Bessel functions is now given by normalizing /i„ 
according to the sum rule (fTO|) . 



J„(x,y) w sgn ( J 



M, 



(30) 



n=M_ 



The reason why we normalize the sum of squares is 
that a summation of only nonnegative terms cannot 
suffer from numerical cancellation. An alternative 
way of normalization would consist in calculating 
a particular value of Jn{x,y), say Jo{x,y), by an- 
other method, like the sum ([5]), or an asymptotic 
expansion [33|. In this case, the approximation to 
Jn{x, y) would be given as 



for all 



Jn{x,y) 



< n < n„ 



Joix,y) 
ho 



hn , 



(31) 



provided Jo(x, y), /iq ^ 0. 



To illustrate some of the intermediate steps of the algo- 
rithm, we show in Fig.[3]the typical behavior of the coef- 
ficients and A^'^ calculated in step 2, and also the 
result after step 3, before normalization of the arrays /„ 
and gn- Concluding the description of our recursive algo- 
rithm, we summarize the different integer indices which 
occur in the problem, which is useful to have in mind 
in the ensuing discussion: n_ and n+ are the negative 
and positive cutoff indices, respectively, and are fixed by 
the values of the arguments x and y through Eq. 
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M_n- n+M+ 




(a) (b) 

FIG. 3: (Color online.) Panel (a) shows the absolute value of the coefficients from Eq. (|27|) and Eq. (I28j . used for the transformed 
recurrence relations (|25|l and (|26p . Here, a starting index of Af_ = —4000 and initial values of = X^j^^ = 1 were used. 

Note that for better visibility, all curves except have been shifted vertically on the logarithmic ordinate axis by multiplication 
with a suitable constant (10^" for ^fj, 10^ for 10~^ for A^, and lO"^'^ for A^). In panel (b), we display the absolute values 
of /„ and Qn, which result after completing step 3 in the algorithm described, i.e. before normalizing. The dash-dotted lines 
indicate the starting indices M±, here M_ = —3300, and M+ = 2350. The cutoff indices n_ = —3000, n+ = 2063 are plotted 
with dashed lines. An example of a suitable "matching index" K = (see step 4 of the algorithm) is drawn by a solid line. 
The initial values used to calculate the curves were /m+ = 10~^''/2, /m^+i ~ 0, /m_|_+2 ~ 10"^°, and gM^ = 0, gM^+i = 10^^. 
Note that in panel (b), no vertical shifting was applied. The inset shows a magnification of the cutoff region for positive n, 
where the diverging behavior of Qn for n > n+ is clearly seen. In both graphs we have x = y — 10^, same as in Fig. (2] 



Af_ and are the negative and positive starting in- 
dices, respectively. For the algorithm to converge, they 
should be chosen such that M_ < n_, and M+ > n_|_. 
The accuracy of the computed approximation to J„(x, y) 
will increase if the distances n_ — M_, M+ — n+ are 
increased (see Sec. ITV)) . K is a matching index, where 
the solutions /„ and </„ computed with different recur- 
rence relations should be matched, and should satisfy 
n_ < K < n+. Finally, nmin and rimax are the indices 
between which numerical values for Jn{x,y) are sought. 
Except for the requirements nmin > ^^max < M^, 

and nmin < ?^inax, they can be arbitrarily chosen. The 
usual requirement is however rimin < rt_, rimax > n^, 
which in that case implies the following inequality chain 
for the different indices involved: 

M_ < n^in <n-<K<n+< n^ax < M+. (32) 



C. Demonstration of numerical stability 

In this subsection we show, using arguments similar 
to those in [36|, that the previously presented algorithm 
is numerically stable. Since the functions J„, Yn, X„ 
and Zn (see Fig. [5]) form a complete set of functions with 
respect to the recurrence relation (llip . we can decompose 
any solution ?/„ to the four-term recurrence relation (|25p 



as 

yn = ai J„ + 02 Yn + fla X„ + 04 Z„. (33) 

The constants ai, 02, 03, 04 can be found from the initial 
conditions. For general i in the range N — 2 < i < N, 
where TV is a general starting index (later we will take 
N = M_), we have 

yi ^ aiJi + a2 Yi + 03 + 04 Zi, (34) 

but we can rewrite yN+i using the four-term recurrence 
in Eq. ^ as 

VN+l {^NyN+ VN-l + yN~2) 

(35) 

= fll Jn+1 + 0.2 Yn+1 + 03 Xat+i + 04 Zat+i, 

for fixed starting integer N. 

If we now for simplicity take the initial value at the up- 
per boundary of the recursion yN+i = 0, by selecting the 
initial values for the coefficients accordingly, then 

we can choose (provided the system (j35p is nonsingular, 
so that a solution exists) three sets of initial conditions 
Hi, ^ ^ t < 3, N — 2 < i < N, so that depending on 
which set is chosen, the constants aj in Eq. ([55)1 are 

flj = Sjt, 1 < j < 3, 1 < t < 3, (36) 
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where Sjt is the Kronecker delta, leading to the solutions 
ul^ with 1 < i < 3. By requiring ([55]) . we have implicitly 
reduced the solution to a linear combination of just two 
solutions, with nonvanishing components of one of J„, 
Xn, Yn on the one hand, and Z„ on the other hand. The 
remaining constant 04 is obtained, for each set, from 

Un+i = = ^uJn+i + S2tYN+i + SsiXn+i + a^ZiM+i , 

(37) 

assuming ^w+i 7^ 0. Because we have reduced the solu- 
tions ul^ to be linear combinations of just two functions, 
we immediately see that the three sets of initial values 
correspond to the three fundamental solutions y^"^'^ to 
the four-term recurrence relation (1251). 



Vn 



Jn 







Zn+1 




Yn+i 


z, 


Zn+1 




Xn+ 


1 r 



(38) 



^N+l 



If now N is taken small enough, N — Af_ < n_, by 
virtue of Eq. (fT9|) . the three fundamental solutions y^^'^ 
turn to the three functions J„, Yn, and X„. We have ba- 
sically eliminated the unwanted solution Z„ by rewriting 
the five-term recurrence into a four-term recurrence 
relation (^5)) . 

In other words, the reduced four-term recurrence rela- 
tion ([25|l . with the coefhcients f,^^''^ evaluated according 
to (j27p in the direction of increasing n from initial val- 
ues ' , < ?i_ with a safety margin, has to a very 
good approximation the three functions J„, y„, and Xn 
as fundamental solutions. This means that a solution /„ 
to the recurrence relation ((25)) . started with initial values 
/;, M+ < I < M+ -I- 2, M+ > n+ with a safety margin, 
and applied in the direction of decreasing n will be almost 
proportional to J„ for n < M^, by the same arguments 
as in subsection IIII Al because after having eliminated 
Zn, the wanted solution J„ is the only one which is sup- 
pressed for n — > 00. This is however only true down to 
the negative cutoff index n_ below which an admixture 
of the other unwanted solution Yn takes over, see Fig. |3l 

Similarly, for the three-term recurrence relation (|26p . 
we can write a generic solution Vn in terms of the three 
fundamental solutions to the four-term recurrence rela- 
tion 111), 



Vn=biJn+b2Xn + hY„. (39) 

Again, there exist two sets Vj,l<s<2, N— l<j< 



iV -|- 1, of initial conditions, with w]^^^ ~ that 

- S,j. (40) 
The two fundamental solutions to Eq. ([^5]) are therefore 

(41) 



- J --Ii^+Iy 



vi^Xn- 



Yn+1 
Xn+i 
Yn+i 



Thus, provided the recurrence for the coefficients of the 
three-term recurrence given in Eq. (|28p is started at suf- 
ficiently small, negative N = M_ < and applied 
in the forward direction, a solution gn to the three-term 
recurrence relation ([26]), started at a large M+ > n+ 
and performed in the direction of decreasing n, will, to a 
good approximation, be proportional to J„ for n < n_. 
Combining the solution /„ to the four-term equation (I25p 
with the solution gn to the three-term equation (|26p at 
the matching index K, where n_ < K < then yields 
a solution proportional to J„ for all n, Timin < < J^max- 
The proportionality constant is found using the sum rule 

Having settled the question of convergence, we com- 
ment briefly on how to numerically calculate the comple- 
mentary solutions Ym X„, and Z„, shown in Fig. [21 We 
assume the most interesting case 1, a; < By. The function 
Xn can be computed by using the original recurrence re- 
lation (jlip in the direction of increasing n, starting at an 
index N < —2y + x, i.e. in region bi. Here Xn quickly 
outgrows the other solutions to leave only the "pure" 
Xn after a few iterations. For Zn, we similarly use the 
original recurrence relation (|lip . but this time in the di- 
rection of decreasing n, and starting at a large positive 
index N > [for the definition of n+, see Eq. ([M)) ]. 
However, in this region Jn grows as fast as Zn, and a 
solution yn calculated this way will be a linear combina- 
tion yn — ai Jn + 0,2 Zn for n > —2y + x, the constants 
ai.2 depending on the initial values. For n < —2y-\-x (in 
region 6i), Z„ grows faster with decreasing n than the 
other fundamental solutions, so that here j/„ — a2Z„. 
Finally, using the four-term relation ([25)1 in the backward 
direction, starting at index N < n+ in region ci , yields a 
solution Xn = aiJn + a'2Yn + a3Xn for — 2y-f x < rt < 7i_|-, 
Xn — aiJn + a2Yn for < n < —2y + x and Xn — a2Yn 
for n < n_. 



IV. DISCUSSION 

A. Accuracy 

It is necessary to investigate how the accuracy of the 
computed approximation /i„ depends on the starting in- 
dices M_, M-|_. To this end, we define the positive "safety 
margin" parameter A through 



M+ ^ n.y 



(42) 



so that specifying A fixes both the upper and the lower 
starting index, and we also define the relative error 



hn - Jn{x,y) 



Jn{x,y) 



(43) 



Yn 



In Fig. [31 we show the relative accuracy that can be 
obtained by the method presented in this paper, as a 
function of A, for different values of the arguments x 
and y, and different index n in the obtained array /i„. 
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We have numerically verified that a performance, simi- 
lar to the one presented in Fig. 01 can be expected even 
close to zeros of J„(a;,y) (that is, for a general index n, 
"mill <n< Timax, whcrc Jn[x, J/) = Or J„(x, y) « 0), al- 
though in this case the estimates remain valid only for the 
absolute instead of the relative error. Specifically, in the 
panels (a)-(c) of Fig. 21 we evaluate the relative error e^c\ 
at n = and take rimin = "max = [see Eq. (1^^ . 
and also the discussion preceding Eq. (|32p ]. which means 
that the recurrence is started at a distance A from the 
cutoff indices. The different curves in the graphs cor- 
respond to the following values of x and y: In (a), we 
have 2y = a; = 10 for curve 1 (blue line), 2y — x = 10^ 
for curve 2 (green line), and 2y = x ~ 10'^ for curve 3 
(red line). For these values of x, y, the index n — Q 
corresponds to the border between the two saddle point 
regions bi and ci. We note that J„(a;, y) cannot be accu- 
rately evaluated in such border regions using the simple 
saddle point approximation [30, 32], but that our method 
works well here. In (b) we have y = IQx — 10 for curve 1 
(blue line), y = lO^cc = 10^ for curve 2 (green line), and 
y = IQ^x = 10^ for curve 3 (red line), demonstrating the 
method for cases where the ratio y/x is large. In (c), we 
have instead a small ratio y/x: x = lOy = 10 for curve 
1 (blue line), x = lO^y = lO'^ for curve 2 (green line), 
and X — lO'^y — 10'^ for curve 3 (red line). Finally, in 
(d) we show the case where Jn{x,y) is evaluated in the 
cutoff region, where for all three curves \ Jn{x,y)\ is of 
order 10"^". Here, we have y = x = 10, n = rimax — 55, 
n,„in = n_ — n -|- n+ = —64 for curve 1 (blue line), y = 
X — 10^, n = TT-niax = 270, rimin = ri_— n + fi 
curve 2 (green line), and y = x = 10^, n — n 
nmin = ri- — n + = —3137 for curve 3 (red line). The 
value nmin has in all cases in graph (d) been chosen so 
that the distance n_ — nmin equals nmax — n^.. Recall that 
the starting indices M± follows by fixing nmin, nmax, and 
A, by Eq. (HH). The black circles in the graphs (a)-(d) 
have been obtained from Eq. (|45p . using approximation 
For the calculations, computer arithmetic with 32 
decimals was used. 

An analytic formula for the relative error can be ob- 
tained by assuming that after normalization, the calcu- 
lated value hn is of the form (writing out the dependence 
of Yn on the arguments x and y explicitly) , 



-364 for 

max ~ 2200, 



Jn{x,y) 



JM^{x,y) 
Ym^ {x, y) 



Yn{x,y), 



(44) 



for starting index M_ < n_, similarly to the case for 
the ordinary Bessel function [see Eq. 1221)] • This is a 
simplified assumption, since the total error in general is 
more complicated, but Eq. ((44l) can nevertheless be used 
to make practical predictions about the dependence of 
Crei on A. Equation (|44p yields for the approximative 
relative error 



^rcl, app 



Jn(x.,y) 



Jn{x,y) 



JM^{x,y) Yn{x,y) 



YM_{x,y) Jn{x,y) 




100 
A 



(a) 



100 200 
A 

(b) 




50 100 
A 

(c) 



50 100 150 
A 

(d) 



FIG. 4: (Color online.) The relative error trci, as defined in 
Eq. (|43|l . as a function of the "safety margin" parameter A 
[see Eq. (|42p ]. In each of the different parameter ranges con- 
sidered, an exponential decrease of the obtained error with 
the safety margin parameter is observed, demonstrating the 
applicability of the recursive method. In (a), we consider pa- 
rameters such that X = 2y, in (b) we have large ratio y/x, 
in (c) small ratio y/x, and in (d) results for the cutoff region 
are presented. Detailed explanation of the parameter regions 
considered is in the text. In all graphs, black circles repre- 
sent the approximation for the relative error obtained from 
Eqs. (|45 l) and (|46l) . 



An approximation for the amplitude of Ym_ [x, y) for 
M_ < n_ can be obtained from the saddle point ex- 
pression IHl) for Jn{x,y), but reversing the sign of the 
real part of the argument of the exponential. If we write 
the saddle point approximation of Yn{x, y) as Yn{x, y) = 
G+{n) -f G-{n), we have 



I Jm_ {x,y)\ 



Jm^ {x,y) 




F+(Af_)- 


f i^_(Af_) 


Ym^ {x, y) 




G+(M_) - 


f G_(Af_) 


« fe-l^° 




+ )-iAf_0+]| 


_^e-|Ro[i/(0 



(45) 



(46) 

where f{9) is defined as after Eq. (fT4ll . and 6± denote the 
two different saddle point solutions from Eq. p^ . with 
n = M_. The last approximation in Eq. ([46l) neglects 
the preexponential factor and the oscillating factor in the 
saddle point approximation (|14p , which is sufficient for an 
order-of-magnitude estimate. The ratio Yn{x,y) / Jn{x,y) 
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in (j45p can be approximated with unity for n in the os- 
cillating region [graph (a), (b), and (c) in Fig. 2], and 
with the simplified saddle point approximation (|46p for 
n in the cutoff region [graph (d) in Fig. 3] . The approx- 
imation (|45p together with ((46|) for the relative error is 
plotted as circles in Fig. [H Clearly the approximate for- 
mula can be used for practical estimates of how far out 
the recurrence should be started if a specific accuracy is 
sought for the array of generalized Bessel functions to be 
computed. Formula (|46p also explains the exponential 
decrease in relative error observed in Fig. 3) 



V. ILLUSTRATIVE CONSIDERATIONS FOR 
THE DIRAC-VOLKOV SOLUTIONS 

In this Section, we consider the Volkov solution, the an- 
alytic solution to the Dirac (or Klein-Gordon) equation 
coupled to an external, plane-wave laser field. We show 
that the generalized Bessel functions can be directly in- 
terpreted as the amplitudes for discrete energy levels of 
a quantum laser-dressed electron, corresponding to the 
absorption or emission of a specific number of laser pho- 
tons. 



B. Comparison with other methods 



A. Physical origin of the recurrence relation 



Here we briefly comment on the performance of the 
presented algorithm as compared to other ways of nu- 
merically evaluating the generalized Bessel function. Let 
us compare to an alternative algorithm based on the eval- 
uation of ordinary Bessel functions using the techniques 
outlined in Sec. IIII Al where we first calculate two arrays 
J2s+n{x), Js{y) of ordinary Bessel functions by Miller's 
algorithm and later calculate the generalized Bessel func- 
tions using Eq. © . Calculation of the arrays of ordinary 
Bessel functions then requires two recurrence runs, and 
to obtain the numerical value J„(x,?/), in addition the 
sum X]^-oo '^2s-i-n(2:)>/s(2/) has to be performed. This 
means that since the generalized Miller's algorithm re- 
quires two recurrence runs only, for calculation of a single 
value J„g{x,y), the two methods demand a comparable 
amount of time. However, the calculation of a single 
generalized Bessel function is not the aim of our consid- 
erations: for the whole array Jn{x,y), nmin < n < Jimax, 
the reduction in computer time due to the elimination of 
the calculation of the sums J27L-oo J2s+n{x)Js{v) leads 
to an order-of-magnitude gain with respect to computa- 
tional resources while the accuracy obtained by the two 
different methods is similar. 

The second method with which to compare is the 
asymptotic expansion by integration through the saddle 
points, as presented in (3^ . For evaluation of a single 
value J„g{x,y), with moderate accuracy demands, the 
saddle-point integration is of course the best method, es- 
pecially for large values of the parameters n, x and y. 
The drawback of this method is the relatively complex 
implementation S^l, and in addition, an increase in the 
accuracy of a saddle-point method typically is a nontriv- 
ial task which involves higher-order expansions of the in- 
tegrand about the saddle point, and this typically leads to 
very complicated analytic expressions for higher orders, 
especially for an integrand with a nontrivial structure as 
in Eq. ([1]) . See however [s^ for a possibly simpler numeri- 
cal method, the "numerical steepest descent method" . In 
any case, if the complete array Jn(a;, y), <n< Timax 
is sought to high accuracy, as it is the case for second- 
order laser-related problems, then our method is neces- 
sarily better, since the time spent on one recursive step 
is very brief. 



There is a direct, physical way to derive the recurrence 
relation satisfied by the generalized Bessel function, in 
the context of relativistic laser-matter interactions. The 
result of this approach defines Jn{x,y) in terms of the 
recurrence relation and a normalization condition, even 
on the level of spinless particles, i.e. on the level of Klein- 
Gordon theory. In this section, we set h = c = 1, denote 
the electron's charge and mass by e = — |e| and m, re- 
spectively, and write dot products between relativistic 
four-vectors a,s u ■ v = u^v^ = u^v^ — u ■ v, for two four- 
vectors and v^^. The space-time coordinate is denoted 
by = (t,x), in order not to cause confusion with the 
argument x of J„ {x,y), and k - z — ujt^ k - x is the phase 
of the laser field. The 4x4 Dirac gamma matrices are 
written as 7^^. 

Let us consider the Klein-Gordon equation 

{idz — eAf — rm? tp{z) = for the interaction of 

a spinless particle of charge e with an external laser field 
of linear polarization A^^{z) = cos{k ■ z), 



9f — 2ie cos(A: ■ z) a ■ dz 



■ cos(2fc • z) — — 



|eV| 



(47) 



^{z) = 0. 



Here a'^ = (0, o) is the polarization vector, and k^^ — 
(w, k) is the propagation wave vector of the laser field, 
with k ■ k — k ■ a — 0. We also introduce the four-vector 
the so-called effective momentum [sl], which fulfills 



2 2 
q = m 



i|eV| 



(48) 



We now insert the Floquet ansatz [54] for the wave func- 
tion 



-\sk-z 



(49) 



where the coefficients i?„ are independent of z^, into 
Eq. (|47p . From this representation, we see that the fac- 
tor e""*''"'^, actually has the same form as a phase factor 
characterizing the absorption of s laser photons from the 
laser field, as we integrate over the Minkowski coordi- 
nate z in the calculation of an S'-matrix element. For 
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negative s, we instead have emission into the laser mode. 
Equation also leads to a relation for the coefficients 



cos((/)) - 2y cos{2(t)) - s] B, e"'^"^ = , 



S — — 00 



(50) 



ea ■ q 
X = , y 



k ■ q 



k ■ z . 



■ q 



Multiplying Eq. ([50|) with e'"'^, and integrating over one 
period, we obtain the recurrence relation pip , if we iden- 
tify 



(51) 



For the wave function constructed from the solution 
to the recurrence relation (jlip to be finite, we must de- 
mand Jn{x,y) to be normalizable. This is expressed by 
the condition (fTU|). Furthermore, using the property ^ 
of Jn{x,y), we can perform the sum over s in with 
the result 



(52) 



which is the form in which the Volkov solution is usually 
presented j53i] . 

For comparison, the solution '^{z) to the Dirac equa- 
tion in presence of a linearly polarized laser field. 



(i7 • 9z — ea • 7 COS0 — m) 'I'(z) = (53) 



reads 



*(z)=e-''- ( ' ;r;^ [J.+i(x,y) + J,_i(x,y)] 



Ak ■ q 

+ Js{x,y)) c-'^^u. 



1 ' 



(54) 



where x, y are defined in Eq. (|50p . and Uq is a Dirac 
bispinor satisfying 



■ q + 2yj ■ k — m) Ug — . 



(55) 



The four-vector p^^ = q'^ + 2yfe^ can be identified as the 
asymptotic momentum of the particle, or the residual 
momentum as the laser field is turned off. 



B. Classical-quantum correspondence of Volkov 
states 



It follows from the expression (j49|l . that a quantum 
Volkov state (we consider the spinless case for simplicity) 
can be regarded as a superposition of an infinite num- 
ber of plane waves with definite, discrete, four-momenta 

-I- nk^^ . The amplitude to find the particle with four- 
momentum + nfc^ is given by Jn{x,y), with x and y 
as in Eq. ([50]) . Therefore, it might seem that the particle 
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FIG. 5: (Color online.) Illustration of the classical-quantum 
correspondence of a laser-dressed electron. Shown to the left 
[panel (a)] with a solid blue line is the classical energy it'' [see 
Eq. (|57|l ] as a function of the phase (f>, for cj — 1 eV, laser 
intensity I = 10^® W/cm^ (corresponding to |ea|/m = 0.1 
where a is the laser polarization four- vector) , initial energy 
p'^ = 2m, and initial angle 6 = 0.54°, with p ■ k — \p\io cos 6. 
These parameters give x — 3.3 x 10"^, and y — —2.7 x 10^ 
for the arguments of J„{x,y). The dashed lines show the 
minimum classical energy Mmim the maximum classical en- 
ergy u^ax) the average energy q'^, and the intermediate en- 
ergy level (local maximum) uf^^, as indicated in the center of 
the figure. In panel (b), we display the quantum mechanical 
amplitude Jn{x, y) of energy level n, which has energy q" + nuj 
[see Eq. (I49p ]. Here positive index n corresponds to absorb- 
ing n number of photons from the laser field, while negative 
n means emitting |n| number of photons into the laser mode. 
The graph is arranged such that n — corresponds to the 
average energy q". The cutoff indices are nicely reproduced 
by the classical maxima and minima. 



can acquire arbitrarily high energy in the field. That this 
is not so, follows from the exponential decay of Jn{x,y) 
beyond the cutoff indices, as discussed in subsection lllBI 
In the following, we show that the cutoff indices can also 
be derived as the lowest and highest energy of a classical 
particle moving in a laser field. To this end, we first recall 
the classical, relativistic equations of motion of a particle 
of charge e and mass to, moving in the laser potential 



a^^ cos ( 



dii^ 6 

— — = — (a^ fc • u - fc'^ a • u) sin(^, (56) 
dr TO 

where is the kinetic momentum, and r is the proper 
time. The solution reads [flU^j assuming initial phase 
= 7r/2, 

u^^ ^ + xk^^ cos - 4 y fc^ cos^ 4>- ea^ cos , (57) 

where is the asymptotic momentum. Note that as u^^ 
is the physical momentum, it is gauge invariant under 
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A^^ + Ak^, where A is an arbitrary function. The 

phase average is exactly the effective momentum, = 
pM _ k^e^a^ /{Ak ■ p) — qf^. In Fig. \5\ we consider the 
energy u° as a function of the phase (p and compare it 
with the discrete energy levels + nuj of the quantum 
wave function. We see that the maximal and minimal 
energy of the classical particle correspond exactly to the 
cutoff indices of the generalized Bessel function. The 
probability for the quantum particle to have an energy 
larger (or smaller) than the classically allowed energy is 
thus exponentially small. Interestingly, the local maxima 
of the classical energy u", labeled uP^^. in Fig.O coincide 
with the transition between the two different saddle point 
regions of J„{x,y). 



VI. CONCLUSIONS 

We have presented a recursive algorithm for numeri- 
cal evaluation of the generalized Bessel function Jn{x, y), 
which is important for laser-physics related problems, 
where the evaluation of large arrays of generalized Bessel 
functions is crucial. In general, we can say that the laser 
parameters fix the arguments x and y of the generalized 
Bessel function Jn{x,y), while the index n characterizes 
the number of exchanged laser photons. 

As evident from Figs. [1] and [H complementary solu- 
tions y„, Xn, and Zn to the recurrence relation (jlip 
satisfied by Jn{x,y) are central to our algorithm. By 
removing the sources of numerical instability, which are 
the exponentially growing complementary solutions, in 
a first recurrence run, we are able to construct a sta- 
ble recursive algorithm, similar to Miller's algorithm for 
the ordinary Bessel function, but suitably enhanced for 
the generalized Bessel function. Numerical stability is 
demonstrated, and the obtainable accuracy is studied nu- 
merically and by an approximate formula (see Sec. lIVp . 
The algorithm is useful especially when a large number of 
generalized Bessel function of different index, but of the 
same argument, is to be generated. As is evident from 
the discussion in Sec. |Vl a fast and accurate calculation 
of generalized Bessel functions leads to a quantitative un- 
derstanding of the quantum-classical correspondence for 
a laser-dressed electron. 
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APPENDIX: INTEGRAL REPRESENTATION OF 
THE COMPLEMENTARY SOLUTIONS 

In this Appendix, we present the expressions for the 
integral representations of the complementary solutions 
Yn{x,y), Zn{x,y), and Xn{x,y) to the recurrence rela- 
tion (jlip . without giving any details about the math- 
ematical considerations which lead to these representa- 
tions. The integrals read 

1 1"°° 

Yn{x,y) = / [cos(n7r/2 - a;cosh6')e"^ 



+ (-l)"e- 



c sinh Oi^—y sinh 29 



/2 



sin(7i6' — x sin 9 + y sin 29)d9, 



(A.l) 



Xn{x,y) 



Zn(x,y) 



1 r°° 

- / sin(n7r/2- a; cosh 6l)e"^-^ "'"'^29(^51 
I" Jo 

1 r/^ 

/ cos(n6' — a;sin6' -I- ?/sin26')d0, 

Jo 

(A.2) 

1 r 

— / sin(n6' — xsin^ + j/sin26')d6' 

Jo 



[(-l)"e 



— a: sinh ^ sinh 6 



^-n9-y sinh 



(A.3) 



Recall that we consider non-zero, positive values of the 
arguments x and ?/, and an arbitrary integer n. By par- 
tial integration, the functions (jA.ip — (jA.3|) verify the re- 
currence relation (fTTjl . The prefactor has been selected 
for each case so that the functions X„(x,y), Yn{x,y), 
and Zn{x,y) have the same amplitude as Jn{x,y) in the 
oscillating region, and this choice also implies that the 
functions Yn{x — > 0,y) and Zn{x 0,y), for even and 
odd n, respectively, can be expressed as Neumann func- 
tions of fractional order. (The latter statement is also 
given here without proof.) A more detailed discussion of 
the mathematical properties of the four functions defined 
by the integral representations ^ , (jA.ip , (|A.2|) and (|A.3|) 
will be given elsewhere. For all considerations reported 
in the current article, the detailed knowledge of the in- 
tegral representations is not necessary; it is sufficient to 
know the recurrence relation (fTTjl that they fulfill. 
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